/*************************************************************************
Purpose:  This program tests various solvers using random numbers.
************************************************************************** 
Input: None.
**************************************************************************
Output: Displays solutions and errors, that's all.
**************************************************************************   
Copyright 2010 Bryan Killett

This file is part of GAIA.

GAIA is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

GAIA is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with GAIA. If not, see <http://www.gnu.org/licenses/>.
*************************************************************************/
//Header file "definitions.hpp" declares many functions, constants, datatypes.
//In order to use the functions declared in definitions.hpp, must link
//with object code that includes functions.cpp.
#include "definitions.hpp"

int main() {

  //These variables are used to display run time.
  time_t start_time,end_time;
  timespan_conversion_s timespan_conversion;//From secs to human units
  
  long long i,j;
  least_squares_s least_squares;
  vector<double> temp_vector_double;//Used to push_back vector of vecs
  vector< vector <double> > results;//First index for solver, second for soln vectors.
  vector<string> titles;
  
  /////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////
  //Define random vector and matrix to be solved via GSL and others.
  start_time = time(NULL);//Record time at start.
  double sigma = 1.0;
  
  //Requires GSL libraries.
  #ifdef GSL_HERE
    //Record the dimension of b (and first dimension of RECTANGULAR A) here.
    least_squares.dimension = 100;
    long long rows = least_squares.dimension;//Num of data pts = dimension (i.e. num of params).
    //long long rows = 10000;//Num of data pts > dimension (i.e. num of params).

    cout<<"Allocating GSL A matrix with dimensions "<<least_squares.dimension<<", "<<rows<<" and GSL b vector with dimension "<<rows<<endl;
    //Initialize "A" matrix- correct dimensions and populated with 0.0.
    least_squares.A = gsl_matrix_alloc(rows, least_squares.dimension);
    gsl_matrix_set_zero(least_squares.A);

    //Initialize "b" vector- correct dimension and populated with 0.0.
    least_squares.b = gsl_vector_alloc(rows);

    gsl_vector_set_zero(least_squares.b);
    const gsl_rng_type * R;
    gsl_rng * r;
    
    gsl_rng_env_setup();
    R = gsl_rng_default;
    r = gsl_rng_alloc(R);
    gsl_rng_set(r,time(0));

    for(i=0;i<rows;i++){
      gsl_vector_set(least_squares.b, i, gsl_ran_gaussian(r,sigma));
      for(j=0;j<least_squares.dimension;j++){
        gsl_matrix_set(least_squares.A, i, j, gsl_ran_gaussian(r,sigma));
      }
    }
      
  #else 
    cout<<"Random matrix generation requires GSL. Search definitions.cpp for GSL_HERE."<<endl;
  #endif

  //Just in case.
  least_squares_s least_squares_backup = least_squares;

  end_time = time(NULL);//Record time at end.
  timespan_conversion.seconds = (long long)(end_time - start_time);
  cout<<"Random matrix generation finished after"<<sec2human(timespan_conversion)<<"."<<endl;

  /////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////
  //Solve via GSL LU decomposition.
  if(1==1){
    start_time = time(NULL);//Record time at start.
    //Just in case.
    least_squares = least_squares_backup;
    results.push_back(temp_vector_double);
    titles.push_back("GSL LU decomp");

    gsl_vector * x = gsl_vector_alloc(least_squares.dimension);
    int gsl_s;
    gsl_permutation * p = gsl_permutation_alloc(least_squares.dimension);
    cout<<"Solving b=Ax using gsl_linalg_LU_decomp() and gsl_linalg_LU_solve().  This may take a while..."<<endl;

    gsl_linalg_LU_decomp(least_squares.A, p, &gsl_s);
    gsl_linalg_LU_solve (least_squares.A, p, least_squares.b, x);

    for(i=0;i<least_squares.dimension;i++){
      results.back().push_back(gsl_vector_get(x,i));
      //cout<<"results["<<results.size()-1<<"]["<<i<<"] = "<<results.back().at(i)<<endl;
    }

    //Free the space taken up by the various matrices and vectors.
    gsl_matrix_free(least_squares.A);
    gsl_vector_free(least_squares.b);
    gsl_vector_free(x);
    gsl_permutation_free(p);

    end_time = time(NULL);//Record time at end.
    timespan_conversion.seconds = (long long)(end_time - start_time);
    cout<<titles.back()<<" finished after"<<sec2human(timespan_conversion)<<"."<<endl;
  }
  /////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////
  //Solve via GSL Householder.
  if(1==1){
    start_time = time(NULL);//Record time at start.
    //Just in case.
    least_squares = least_squares_backup;
    results.push_back(temp_vector_double);
    titles.push_back("GSL HH");

    gsl_vector * x = gsl_vector_alloc(least_squares.dimension);
    cout<<"Solving b=Ax using gsl_linalg_HH_solve().  This may take a while..."<<endl;

    gsl_linalg_HH_solve(least_squares.A, least_squares.b, x);

    for(i=0;i<least_squares.dimension;i++){
      results.back().push_back(gsl_vector_get(x,i));
      //cout<<"results["<<results.size()-1<<"]["<<i<<"] = "<<results.back().at(i)<<endl;
    }

    //Free the space taken up by the various matrices and vectors.
    gsl_matrix_free(least_squares.A);
    gsl_vector_free(least_squares.b);
    gsl_vector_free(x);

    end_time = time(NULL);//Record time at end.
    timespan_conversion.seconds = (long long)(end_time - start_time);
    cout<<titles.back()<<" finished after"<<sec2human(timespan_conversion)<<"."<<endl;
  }
  /////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////
  //Solve via GSL QR decomposition.
  if(1==2){
    start_time = time(NULL);//Record time at start.
    //Just in case.
    least_squares = least_squares_backup;
    results.push_back(temp_vector_double);
    titles.push_back("GSL QR decomp");

    //residual is "e" in b=Ax+e.
    gsl_vector * residual = gsl_vector_alloc(rows);
    gsl_vector * x = gsl_vector_alloc(least_squares.dimension);
    gsl_vector_set_zero(residual);
    //The vector tau contains some kind of internal QR decomp information.
    gsl_vector * tau = gsl_vector_alloc(least_squares.dimension);
    gsl_vector_set_zero(tau);
    cout<<"Solving b=Ax using gsl_linalg_QR_decomp() and gsl_linalg_QR_lssolve().  This may take a while..."<<endl;

    gsl_linalg_QR_decomp(least_squares.A, tau);
    gsl_linalg_QR_lssolve(least_squares.A, tau, least_squares.b, x, residual);

    for(i=0;i<least_squares.dimension;i++){
      results.back().push_back(gsl_vector_get(x,i));
      //cout<<"results["<<results.size()-1<<"]["<<i<<"] = "<<results.back().at(i)<<endl;
    }

    //Free the space taken up by the various matrices and vectors.
    gsl_matrix_free(least_squares.A);
    gsl_vector_free(least_squares.b);
    gsl_vector_free(x);
    gsl_vector_free(tau);
    gsl_vector_free(residual);

    end_time = time(NULL);//Record time at end.
    timespan_conversion.seconds = (long long)(end_time - start_time);
    cout<<titles.back()<<" finished after"<<sec2human(timespan_conversion)<<"."<<endl;
  }
  /////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////
  //Solve via another method.
  if(1==2){
    start_time = time(NULL);//Record time at start.
    //Just in case.
    least_squares = least_squares_backup;
    results.push_back(temp_vector_double);
    titles.push_back("another method");

    end_time = time(NULL);//Record time at end.
    timespan_conversion.seconds = (long long)(end_time - start_time);
    cout<<titles.back()<<" finished after"<<sec2human(timespan_conversion)<<"."<<endl;
  }
  /////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////
  //Compare answers- they should be identical!
  double error, max_error = -1000000;
  for(i=0;i<least_squares.dimension;i++){
    double max=-1000000,min=-max;
    for(j=0;j<(long long)results.size();j++){
      if(results.at(j).at(i) > max) max = results.at(j).at(i);
      if(results.at(j).at(i) < min) min = results.at(j).at(i);
    }
    error = fabs((max-min)/min);
    if(error > max_error) max_error = error;
    //cout<<"Param #"<<i+1<<" has max error "<<error*100.0<<"%."<<endl;
  }
  cout<<"Max error: "<<max_error*100.0<<"%."<<endl;

  return 0;
}
